The following script and analyses were written by Colleen Bove and James Umbanhowar. This version of the script was last committed to GitHub on 22 Febuary 2020

Results

In total, eleven studies met the standards of this meta-analysis, including the responses of thirteen Caribbean coral species collected from five different countries across the wider Caribbean (Figure 1 and Supplementary Table S1). Of the studies selected, only four performed fully factorial ocean acidification and warming experiments, and one performed two independent acidification and warming experiments. The most studied coral species from the Caribbean region were Porites astreoides (5 studies), Acropora cervicornis (4 studies), and Siderastrea siderea (4 studies). Finally, the Florida Keys and Belize were the two most-studied regions within the wider Caribbean fitting the criteria of this meta-analysis.

3.1 Overall calcification response

Meta-analysis of the dataset revealed that calcification rates of Caribbean corals were reduced by ocean warming but not ocean acidification (Figure 2 and Supplementary Figures S2). However, the 95% confidence interval of the combination of ocean warming and acidification overlapped zero, indicating no statistically clear trend towards synergistic or antagonistic effects of these treatments (Figure 2 and Supplementary Table S2).

Figure 1 Coral collection sites of all included studies with experimental study represented by colour and treatments represented by shape: acidification only (circle), warming only (triangle), and the combination of acidification and warming (square). The lower left insert displays close-up of the Belize collection sites and the upper right insert displays the Florida Keys collection sites.

# calculate ES and variance for ocean acidification
ES.oa <- escalc(measure = "SMDH", m1i = cm1, m2i = am3, sd1i = s1, sd2i = s3, 
    n1i = n1, n2i = n3, data = df)

# calculate ES and variance ES for ocean warming
ES.ow <- escalc(measure = "SMDH", m1i = cm1, m2i = wm4, sd1i = s1, sd2i = s4, 
    n1i = n1, n2i = n4, data = df)
ES.ow <- ES.ow %>% rename(yi.w = yi, vi.w = vi)

# create DF of calculated variables of OA
acid <- data.frame(cm1 = ES.oa$cm1, yi.a = ES.oa$yi, vi.a = ES.oa$vi)

# merge together warming and acid DF
ES.full <- merge(ES.ow, acid, by = "cm1")
ES.full <- ES.full[order(ES.full$cite), ]
ES.full$yi.w <- (ES.full$yi.w) * -1
ES.full$yi.a <- (ES.full$yi.a) * -1


### Standard deviation

# calculate sd for all
ES.full$sd1 <- ES.full$s1 * sqrt(ES.full$n1)
ES.full$sd3 <- ES.full$s3 * sqrt(ES.full$n3)
ES.full$sd4 <- ES.full$s4 * sqrt(ES.full$n4)
ES.full$sd6 <- ES.full$s6 * sqrt(ES.full$n6)

# calculate S for all treatments
ES.full$Sa <- sqrt((ES.full$sd1^2/ES.full$n1) + (ES.full$sd3^2/ES.full$n3))
ES.full$Sw <- sqrt((ES.full$sd1^2/ES.full$n1) + (ES.full$sd4^2/ES.full$n4))
ES.full$Si <- sqrt((ES.full$sd1^2/ES.full$n1) + (ES.full$sd6^2/ES.full$n6))


### Effect Size of Interaction term

# calculate ES using Harvey et al equations for acidification, warming, and
# interaction studies
ES.full$lnRRa <- ((ES.full$am3 - ES.full$cm1)/ES.full$Sa)
ES.full$lnRRt <- ((ES.full$wm4 - ES.full$cm1)/ES.full$Sw)
ES.full$yi.i <- ((ES.full$m6 - ES.full$wm4) - (ES.full$am3 - ES.full$cm1))/(2 * 
    ES.full$Si)

# calculate V for interaction studies ONLY
ES.full$vi.i <- (1/ES.full$n1) + (1/ES.full$n3) + (1/ES.full$n4) + (1/ES.full$n6) + 
    ((ES.full$yi.i^2)/(2 * (ES.full$n1 + ES.full$n3 + ES.full$n4 + ES.full$n6)))


write.csv(ES.full, file = "~/Bove2020_MetaAnalysis_FrontMarSci/Files/ES.full.csv")

# read in the saved CSV and then removed any unneeded columns
df<-read.csv("~/Bove2020_MetaAnalysis_FrontMarSci/Files/ES.full.csv")
df<-df[,-c(2, 12:13, 16:18, 20:22, 24:26, 41:49)] #remove any unneeded columns

### data manipulation

#subset df for variance (v) values only as new df
df2b<-subset(df, select=c(X,vi.i,vi.a,vi.w))

#subset df to remove only variance (v) values as another new df
df.m <-subset(df, select=-c(vi.i,vi.a,vi.w))

# this creates df with treat and V as column headers in long format
df2b<-gather(df2b, treat, v,vi.i:vi.w)

#substitute y for v for merge later
df2b$treat<-gsub('v','y', df2b$treat)

#creates another long format df with treat and Y as headers (of Y values)
df.m <-gather(df.m, treat, y,yi.w:yi.i)
df.m $treat<-as.factor(df.m $treat)

#merge two df into one
df2<-merge(df.m, df2b, by=c("X", "treat"))

#remove any NAs
#df2<-na.omit(df2)
df2$study <- paste(df2$author, df2$year, sep = "_")
df2$X<-NULL


## centering treatments

df2$mag.p <- df2$a.pco2 - df2$c.pco2
df2$mag.t <- df2$w.temp - df2$c.temp

df2$mag.cent.p <- df2$mag.p - 127 # make a centered around 0 column
df2$mag.cent.t <- df2$mag.t - 1.0 # make a centered around 0 column

# subset data by treatment
OW<-subset(df2, treat=="yi.w")
OA<-subset(df2, treat=="yi.a")
IN<-subset(df2, treat=="yi.i")

# fit model to each treatment
ow <- rma(y, v, data=OW)
oa <- rma(y, v, data=OA)
both <- rma(y, v, data=IN)

# wide to long format df
df2 <- gather(data = df2, key = "treat.mag", value = "mag", mag.p:mag.t)

# remove NA values from df
completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}

df2 <- completeFun(df2, "y")
df2 <- completeFun(df2, "mag")

global.mod <- rma.mv(y,v, mods = ~ treat , random = list(~ (treat|study), ~ (treat|spec)), data = df2) 

preds <- as.data.frame(predict(global.mod, newdata = df2))[,c(1,3:4)] # pull the mean and CI with the predict function
global.out <- cbind(df2, preds) # add columns to existing df for global output

#global.mod

global.out$treat <- factor(global.out$treat, levels = c("yi.a","yi.w","yi.i"))
global.out$wt <- 1/(global.out$v)

global_CI <- 
  ggplot()+
  theme(panel.grid.major=element_blank(), legend.position="none", panel.grid.minor=element_blank(), panel.background=element_blank())+
  theme(axis.line=element_line(color="black"))+
  theme(legend.key=element_blank())+
  ylab("Mean effect size (SMD)")+
  xlab("")+
  #ggtitle("Figure 2. Global standardized mean difference by treatment")+
  scale_x_discrete(labels=c("yi.a" = "acidification", "yi.w" = "warming", "yi.i" = "combination"))+
  geom_hline(yintercept=0, linetype=2, colour="black")+
  geom_point(aes(x = treat, y = y, size = wt),  colour="grey", data = global.out, alpha = 0.5, position=position_jitter(width=0.1))+
  #geom_errorbar(aes(ymin=lb, ymax=ub, width=0.1), size=0.6)+
  geom_linerange(data = global.out, aes(x = treat, ymin=ci.lb, ymax=ci.ub), size=0.8)+
  #ylim(-1.6,0.6)+
  geom_point(data = global.out, aes(x = treat, y = pred), size=10, shape=95)

div(ggplotly(global_CI), align = "center")

ggsave("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Manuscript/Figure2_GlobalCI.pdf", width = 5, height = 3.5)

Figure 2 Mean effect (standard mean difference) and 95% confidence interval of ocean acidification, warming, and the combination of acidification and warming on calcification rate for all studies in the meta-analysis. Grey circles indicate the effect size of each individual study and the size of each circle represents the weight of each study (1/SE). Clear statistical evidence of a treatment effect is identified when the 95% confidence interval does not overlap zero.

3.2 Calcification response of Florida Keys versus Belize corals

Corals from Belize only exhibited clearly reduced calcification rates under ocean warming (Figure 3A and Supplementary Table 3), while acidification, warming, and the combination of both stressors did not clearly alter experimental calcification rates of corals from the Florida Keys (Figure 3B and Supplementary Table S3). Further, the resulting Q_E suggests there is significant between-study variation (Supplementary Tables S3).

df2a<-subset(df2, region!="Curacao")
df2a<-subset(df2a, region!="Little Cayman")
df2a<-subset(df2a, region!="Mexico")

region.mod <- rma.mv(y,v, mods = ~ treat * region, random = list(~ (treat|study), ~ (treat|spec)), data = df2a)

preds2 <- as.data.frame(predict(region.mod, newdata = df2a))[,c(1,3:4)] # pull the mean and CI with the predict function
region.out <- cbind(df2a, preds2) # add columns to existing df for global output

#region.mod

region.out$treat <- factor(region.out$treat, levels = c("yi.a","yi.w","yi.i"))
region.out$wt <- 1/(region.out$v)

region_CI <- ggplot()+
  theme(panel.grid.major=element_blank(), legend.position="none", panel.grid.minor=element_blank(), panel.background=element_blank())+
  theme(axis.line=element_line(color="black"))+
  theme(legend.key=element_blank())+
  ylab("Mean effect size (SMD)")+
  xlab("")+
  scale_x_discrete(labels=c("yi.a" = "acidification", "yi.w" = "warming", "yi.i" = "combination"))+
  geom_hline(yintercept=0, linetype=2, colour="black")+
  geom_point(aes(x = treat, y = y, size = wt),  colour="grey", data = region.out, alpha = 0.5, position=position_jitter(width=0.1))+
  #ggtitle("Figure 3.Standardized mean difference by treatment for Florida Keys and Southern Belize")+
  #geom_errorbar(aes(ymin=lb, ymax=ub, width=0.1), size=0.6)+
  geom_linerange(data = region.out, aes(x = treat, ymin=ci.lb, ymax=ci.ub), size=0.8)+
  geom_point(data = region.out, aes(x = treat, y = pred), size=10, shape=95)+
  facet_wrap(~region)  

div(ggplotly(region_CI), align = "center")

ggsave("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Manuscript/Figure3_RegionalCI.pdf", width = 8, height = 3.5)

Figure 3 Mean effect (standard mean difference) and 95% confidence interval of ocean acidification, warming, and the combination of acidification and warming on calcification rate for (A) Florida Keys corals and (B) Belize corals. Grey circles indicate the effect size of each individual study and the size of each circle represents the weight of each study (1/SE). Clear statistical evidence of a treatment effect is identified when the 95% confidence interval does not overlap zero.

3.3 Temperature and aragnoite saturation state impacts on calcification rates across studies

Secondary analysis of mean calcification rates (mg cm⁻² day⁻¹) against treatment temperature across all Florida and Belize studies revealed a parabolic response to temperature (Figure 4 and Supplementary Tables S4, S5). Similarly, mean calcification rates across Ω_Arag resulted in a nonlinear response to acidification (Figure 5 and Supplementary Tables S6, S7). Both nonlinear trends in response to temperature and Ω_Arag were a result of treatment rather than region (Supplementary Tables S6, S7), suggesting regional differences identified in the meta-analysis were due to experimental designs employed to represent current regional environmental differences.

calc.w<-calc[,-c(15:18, 23:27)] #remove columns not needed for warming analysis

# remove and NA columns from warming column
calc.w <- completeFun(calc.w, "wm4")

# remove mexico studies
calc.w<-subset(calc.w, region != "Mexico")

# subset calc for calcication rates
calc.b<-subset(calc.w, select=c(num, cm1, wm4))
# create df with treat and rate as column headers in long format
calc.b<-gather(calc.b, treat, rate, cm1:wm4)
# rename for useful treatment names
calc.b$treat<-gsub('cm1','control', calc.b$treat)
calc.b$treat<-gsub('wm4','warm', calc.b$treat)

#subset calc for SE
calc.c<-subset(calc.w, select=c(num, s1, s4))
# create df with treat and SE as column headers in long format
calc.c<-gather(calc.c, treat, SE, s1:s4)
# rename for useful treatment names
calc.c$treat<-gsub('s1','control', calc.c$treat)
calc.c$treat<-gsub('s4','warm', calc.c$treat)

#subset calc to remove only columns from above
calc.d <-subset(calc.w, select=-c(s1 ,s4 ,cm1, wm4, c.pco2, n1, n4))
# this creates df with treat and V as column headers in long format
calc.d<-gather(calc.d, treat, temp, c.temp:w.temp)
# rename for useful treatment names
calc.d$treat<-gsub('c.temp','control', calc.d$treat)
calc.d$treat<-gsub('w.temp','warm', calc.d$treat)

#merge dfs into one, final warming
calc.m<-merge(calc.d, calc.b, by=c("num", "treat"))
calc.warm<-merge(calc.m, calc.c, by=c("num", "treat"))

# create column for weight
calc.warm$wt <- 1 / calc.warm$SE

# create study column
calc.warm$study <- paste(calc.warm$author, calc.warm$year, sep = "_")

# center temps
calc.warm$cent.temp <- calc.warm$temp - mean(calc.warm$temp)
calc.warm$cent.pco2 <- rep("NA", 60)

calc.warm$scale.temp <- scale(calc.warm$temp, center = TRUE, scale = TRUE)

warm.rate.region.mod <- lmer(rate ~ scale.temp + region + days + (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F) #sing fit
warm.rate.region2.mod <- lmer(rate ~ scale.temp + region + I(scale.temp^2) + days + (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F) 
warm.rate2.mod <- lmer(rate ~ scale.temp + I(scale.temp^2) + days + (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F) 
warm.rate.mod <- lmer(rate ~ scale.temp + days+ (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F) 
warm.region.mod <- lmer(rate ~ region + days + (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F) #sing fit


#summary(warm.rate.region2.mod) # view summary output of model
#confint(warm.rate.region2.mod, method="boot")

warm.pred <- data.frame(temp=seq(26,32.01,length.out=100), region=rep(c("Belize", "Florida Keys"), length.out = 100), pred=predict(warm.rate.region2.mod, newdata = data.frame(scale.temp=seq(min(calc.warm$scale.temp),max(calc.warm$scale.temp),length.out=100), days=seq(min(calc.warm$days),max(calc.warm$days),length.out=100), region=rep(c("Belize", "Florida Keys"), length.out = 100)), re.form=NA))

global_OW_calc <-
ggplot(calc.warm, aes(x=temp, y=rate))+
  geom_point(aes(size=1/SE, shape=region, fill=region), colour="black")+
    theme(panel.grid.major=element_blank(), legend.position="right", panel.grid.minor=element_blank(), panel.background=element_blank())+
  theme(axis.line=element_line(color="black"))+
  theme(legend.key=element_blank())+
  scale_shape_manual("Region", values = c(21,22,24))+
  scale_fill_manual("Region", values = c("#5ab4ac", "#d8b365", "#f5f5f5"))+
  scale_color_manual("Region", values = c("#5ab4ac", "#d8b365", "#f5f5f5"))+
  geom_line(data = warm.pred, aes(y = pred, group = region, colour = region), size = 1)+
  ylab(bquote('Calcification rate (mg' ~cm^-2* ~day^-1*')'))+
  xlab(bquote('Temperature ('*degree*C*')'))
global_OW_calc

ggsave("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Manuscript/Figure4_CalcWarm.pdf", width = 7, height = 4)

Figure 4. Mean calcification rate (mg cm⁻² day⁻¹) of corals from each warming study by treatment temperature (°C) with the linear mixed effects model quadratic fit. Shape and colour of each point denotes study region (blue circle = Belize; brown square = Florida) and size of shape represents the weight (1/SE) of study.

calc.a<-calc[,-c(19:27)] #remove columns not needed for acid analysis
calc.a$p.mag <- calc.a$c.pco2 - calc.a$a.pco2 # pCO2 magnitude
calc.a <- completeFun(calc.a, "am3") # subset only for acidification studies 
calc.a <- subset(calc.a, region != "Curacao")


#subset calc for calcication rates
calc.b<-subset(calc.a, select=c(num, cm1, am3))
# create df with treat and rate as column headers in long format
calc.b<-gather(calc.b, treat, rate, cm1:am3)
# rename for useful treatment names
calc.b$treat<-gsub('cm1','control', calc.b$treat)
calc.b$treat<-gsub('am3','acid', calc.b$treat)

#subset calc for SE
calc.c<-subset(calc.a, select=c(num, s1, s3))
# create df with treat and SE as column headers in long format
calc.c<-gather(calc.c, treat, SE, s1:s3)
# rename for useful treatment names
calc.c$treat<-gsub('s1','control', calc.c$treat)
calc.c$treat<-gsub('s3','acid', calc.c$treat)

#subset calc to remove only columns from above
calc.d <-subset(calc.a, select=-c(s1 ,s3 ,cm1, am3, c.temp, n1, n3))
# this creates df with treat and V as column headers in long format
calc.d<-gather(calc.d, treat, acid, c.pco2:a.pco2)
# rename for useful treatment names
calc.d$treat<-gsub('c.pco2','control', calc.d$treat)
calc.d$treat<-gsub('a.pco2','acid', calc.d$treat)

#merge dfs into one, final acid
calc.m<-merge(calc.d, calc.b, by=c("num", "treat"))
calc.arag<-merge(calc.m, calc.c, by=c("num", "treat"))

# create column for weight
calc.arag$wt <- 1 / calc.arag$SE

# create study column
calc.arag$study <- paste(calc.arag$author, calc.arag$year, sep = "_")

# center temps
calc.arag$cent.pco2 <- calc.arag$acid - mean(calc.arag$acid)
#calc.arag$cent.pco2 <- rep("NA", 66)

calc.arag$scale.arag <- scale(calc.arag$acid, center = TRUE, scale = TRUE)

acid.rate.region.mod <- lmer(rate ~ scale.arag + region + days + (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F) #sing fit
acid.rate.region2.mod <- lmer(rate ~ scale.arag + region + I(scale.arag^2) + days + (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F) 
acid.rate2.mod <- lmer(rate ~ scale.arag + I(scale.arag^2) + days + (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F) 
acid.rate.mod <- lmer(rate ~ scale.arag + days+ (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F) 
acid.region.mod <- lmer(rate ~ region + days + (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F) #sing fit

#summary(acid.rate.region2.mod) # continuing with this one for simplicity
#confint(acid.rate.region2.mod, method="boot")

acid.pred <- data.frame(acid=seq(0.7,5.8,length.out=100), region=rep(c("Belize", "Florida Keys"), length.out = 100),pred=predict(acid.rate.region2.mod, newdata = data.frame(scale.arag=seq(min(calc.arag$scale.arag),max(calc.arag$scale.arag),length.out=100), days=seq(min(calc.arag$days),max(calc.arag$days), length.out=100), region=rep(c("Belize", "Florida Keys"), length.out=100)),re.form=NA))

global_OA_calc <-
  ggplot(calc.arag, aes(x=acid, y=rate))+
  geom_point(aes(size=1/SE, shape=region, fill=region), colour="black")+
  theme(panel.grid.major=element_blank(), legend.position="right", panel.grid.minor=element_blank(), panel.background=element_blank())+
  theme(axis.line=element_line(color="black"))+
  theme(legend.key=element_blank())+
  scale_shape_manual("Region", values = c(21,22,24))+
  scale_fill_manual("Region", values = c("#5ab4ac", "#d8b365", "#f5f5f5"))+
  scale_color_manual("Region", values = c("#5ab4ac", "#d8b365", "#f5f5f5"))+
  geom_line(data = acid.pred, aes(y = pred, group = region, colour = region), size = 1)+
  ylab(bquote('Calcification rate (mg' ~cm^-2* ~day^-1*')'))+
  xlab(bquote('Aragonite Saturation State'))
global_OA_calc

ggsave("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Manuscript/Figure5_CalcAcid.pdf", width = 7, height = 4)

Figure 5. Mean calcification rate (mg cm⁻² day⁻¹) of corals from each acidification study by treatment aragonite saturation state (Ω_Arag) with the linear mixed effects model quadratic fit. Shape and colour of each point denotes study region (blue circle = Belize; brown square = Florida Keys) and size of shape represents the weight (1/SE) of study.

3.4 Experimental design impacts on coral calcification rate in studies

Quantification of experimental design parameters within warming studies identified that magnitude of treatment, irradiance, seawater type used, and feeding frequency all clearly impacted calcification rates (Table 1). Specifically, studies that utilized natural seawater and those with a larger difference between the control and treatment temperatures within a study exhibited higher effect sizes, suggesting a less negative effect of treatment. Studies that employed higher irradiance levels in their systems demonstrated more negative effects of treatment light level. Finally, studies that reported feeding their corals twice a week were less impacted by warming treatment than those feeding three times a week, however, studies with no data on feeding were the least affected by treatment. Duration of experiment was deemed redundant in the model and was thus dropped.

Within the acidification studies, irradiance, seawater type used, feeding frequency, duration, and the interaction of duration with treatment magnitude impacted effect sizes, while magnitude alone was not clearly different (Table 1). Studies using natural seawater, employing higher irradiance levels, and those with longer durations resulted in great effect sizes, suggesting they lessened the effects of acidification treatment on calcification responses. Similar to warming studies, acidification studies in which feeding corals was conducted twice a week exhibited less negative responses to treatment than those feeding three times a week, with studies reporting no feeding data exhibiting the least negative responses to treatment. Finally, coral calcification responses were less impacted by acidification in studies with longer duration of exposure and a greater pCO2 change.

Table 1. Effect size estimate and 95% confidence intervals of experimental design parameters calculated for all wamring and acidification experiments.

full.param <- rbind(w.global.out[,-4], a.global.out[,-4])

param.tab <- kable(full.param, digits = 3) %>% 
  kable_styling(font_size = 14, full_width = F) %>% 
  pack_rows("Warming Experiments", 1, 6) %>%
  pack_rows("Acidification Experiments", 7, 13)
param.tab

	Estimate	Lower 95% CI	Upper 95% CI	col
Warming Experiments
Intercept	-29.392	-39.564	-19.221	negative
Magnitude	10.341	5.580	15.101	positive
Irradiance	-0.237	-0.336	-0.137	negative
Seawater	29.767	19.719	39.815	positive
Feeding (3x)	-3.544	-4.967	-2.120	negative
Feeding (N.D.)	87.160	49.685	124.635	positive
Acidification Experiments
Intercept	-2.033	-2.785	-1.281	negative
Magnitude	0.000	-0.002	0.003	positive
Irradiance	0.009	0.004	0.014	positive
Seawater	2.392	1.144	3.640	positive
Feeding (3x)	-15.042	-22.607	-7.477	negative
Feeding (N.D.)	3.611	0.280	6.942	positive
Duration	0.206	0.090	0.323	positive
Magnitude x Duration	-0.001	-0.001	0.000	negative

Supplemental Materials

Supplemental Tables

Supplementary Table S1. List of included studies in meta-analysis with experimental design parameters. Treatment temperature and pCO₂ reported are for treatment levels used in this analysis. Missing data from studies is denoted by N.D.

table1 <- kable(meta_tab2, booktabs = T) %>%
  kable_styling(font_size = 10, full_width = F) %>% 
  column_spec(5, italic = T)
table1

	Study	Location	Collection.Sites	Species	Treatments	Experiment.Duration	Feeding	Seawater	Control.Temperature	Temperature	Control.pCO2	pCO2	irradiance
1	Albright et al. 2011	Florida	Little Grecian, Key Largo	Porites astreoides	pCO2	49 d	NA	Natural seawater	28.0	NA	380	800.0	NA
2	Bedwell-Ivers et al. 2017	Little Cayman	Little Cayman Research Center	Acropora cervicornis Porites divaricata	pCO2	28 d	NA	Natural seawater	30.0	NA	487	791.0	800
3	Bove et al. 2019	Belize	Inshore, Offshore (Sapadilla Cayes)	Porites astreoides, Pseudodiploria strigosa, Siderastrea siderea, Undaria tenuifolia	Warming, pCO2, combination	93 d	Frozen and newly hatched Artemia every other day	Filtered non-native seawater	28.0	31.0	405	701.0	300
4	Castillo et al. 2014	Belize	Nearshore, Backreef, Forereef (Sapadilla Cayes)	Siderastrea siderea	Warming, pCO2	95 d	Frozen Artemia every other day	Artificial seawater	28.0	32.0	472	704.0	250
5	Enochs et al. 2014	Florida	Dade County	Acropora cervicornis	pCO2	42 d	NA	Natural seawater	28.0	NA	513	870.9	178
6	Grottoli et al. 2014	Mexico	Puerto Morelos	Porites divaricata, Porites astreoides, Orbicella faveolata	Warming	15 d	Fed	Natural seawater	30.6	31.6	NA	NA	600
7	Horvath et al. 2016	Belize	Nearshore, Backreef, Forereef (Sapadilla Cayes)	Siderastrea siderea	Warming, pCO2, combination	60 d	Frozen Artemia twice weekly	Artificial seawater	28.0	31.0	426	890.0	250
8	Jury et al. 2010	Curacao	Seaquarium	Madracis auretenra	pCO2	2 h	Artemia nauplii twice weekly	Natural seawater	28.0	NA	391	1480.0	200
9	Okazaki et al. 2017	Florida	Key Biscayne, Biscayne National Park, Key West	Acropora cervicornis, Agaricia agaricites, Dichocoenia stokesii, Montastraea cavernosa, Orbicella faveolata, Porites astreoides, Porites divaricata, Pseudodiploria strigosa, Siderastrea siderea, Siderastrea radians, Solenastrea hyades	Warming, pCO2, combination	6 w	Live rotifers and larval twice weekly	Natural seawater	27.0	30.3	400	1340.0	327
11	Towle et al. 2015	Florida	South Florida	Acropora cervicornis	Warming, pCO2, combination	8 w	Dried zooplankton twice weekly	Natural seawater	26.0	30.0	390	800.0	353
12	Kenkel et al. 2013	Florida	Inshore and Offshore Sugarloaf Key	Porites astreoides	Warming	43 d	NA	Natural seawater	27.2	30.9	NA	NA	NA

#save_kable(table1, "Table1_Studies_16Jan20.pdf")

Supplementary Table S2. Global meta-analysis mixed effects model (function rma.mv) output by treatment only, with random effects of study and species, used in Figure 2. The test for residual heterogeneity and significance is represented by Q_E (P-value) and the test of moderators is Q_M (P-value).

tidy.rma <- function(x, ...) {
  with(
    x,
    data.frame(
      estimate = b,
      std.error = se,
      z.value = zval,
      p.value = pval,
      conf.low = ci.lb,
      conf.high = ci.ub
    )
  )
}
tidy.rma.Q <- function(x, ...) {
  with(
    x,
    data.frame(
      QE = QE,
      p.value = QEp,
      QM = QM,
       p.value = QMp
    )
  )
}

global.Q <- tidy.rma.Q(global.mod)
global.Q <- data.frame(t(global.Q))
row.names(global.Q) <- c("QE", "p value", "QM", " p value")
colnames(global.Q) <- "estimate"
global.Q$std.error <- rep("", 4)
global.Q$z.value <- rep("", 4)
global.Q$p.value <- rep("", 4)
global.Q$conf.low <- rep("", 4)
global.Q$conf.high <- rep("", 4)


global.mod.tab <- tidy.rma(global.mod)
#row.names(global.mod.tab) <- c("acidification", "combination", "warming")
global.mod.tab$std.error <- round(global.mod.tab$std.error, 3)
global.mod.tab$z.value <- round(global.mod.tab$z.value, 3)
global.mod.tab$p.value <- round(global.mod.tab$p.value, 3)
global.mod.tab$conf.low <- round(global.mod.tab$conf.low, 3)
global.mod.tab$conf.high <- round(global.mod.tab$conf.high, 3)
global.mod.tab <- rbind(global.mod.tab, global.Q)

global.mod.tab <- global.mod.tab[,1:4]

table2 <- kable(global.mod.tab, digits = 2) %>%
  kable_styling(font_size = 14, full_width = F) %>% 
  pack_rows("Model Results", 1, 3) %>%
  pack_rows("Variance Components", 4, 7)
table2

	estimate	std.error	z.value	p.value
Model Results
intrcpt	-1.20	0.776	-1.545	0.122
treatyi.i	1.19	1.41	0.843	0.399
treatyi.w	-1.07	1.19	-0.898	0.369
Variance Components
QE	823.95
p value	0.00
QM	2.40
p value	0.30

#save_kable(table2, "Table2_GlobalMod_16Jan20.pdf")

Supplementary Table S3. Meta-analysis mixed effects model (function rma.mv) output of treatment by region (Belize versus Florida Keys), with random effects of study and species, used in Figure 3. The test for residual heterogeneity and significance is represented by Q_E (P-value) and the test of moderators is Q_M (P-value).

regional.Q <- tidy.rma.Q(region.mod)
regional.Q <- data.frame(t(regional.Q))
row.names(regional.Q) <- c("QE", "p value", "QM", " p value")
colnames(regional.Q) <- "estimate"
regional.Q$std.error <- rep("", 4)
regional.Q$z.value <- rep("", 4)
regional.Q$p.value <- rep("", 4)
regional.Q$conf.low <- rep("", 4)
regional.Q$conf.high <- rep("", 4)


region.mod.tab <- tidy.rma(region.mod)
#row.names(region.mod.tab) <- c("acidification", "combination", "warming","acidification ", "combination ", "warming ")
region.mod.tab$std.error <- round(region.mod.tab$std.error, 3)
region.mod.tab$z.value <- round(region.mod.tab$z.value, 3)
region.mod.tab$p.value <- round(region.mod.tab$p.value, 3)
region.mod.tab$conf.low <- round(region.mod.tab$conf.low, 3)
region.mod.tab$conf.high <- round(region.mod.tab$conf.high, 3)
region.mod.tab <- rbind(region.mod.tab, regional.Q)

region.mod.tab <- region.mod.tab[,1:4]

table3 <- kable(region.mod.tab, digits = 4, booktabs = T) %>%
  kable_styling(font_size = 14, full_width = F) %>% 
  pack_rows("Belize", 1, 3) %>%
  pack_rows("Florida", 4, 6) %>%
  pack_rows("Variance Components", 7, 10)
table3

	estimate	std.error	z.value	p.value
Belize
intrcpt	-0.9288	1.324	-0.701	0.483
treatyi.i	1.5625	2.24	0.697	0.486
treatyi.w	-3.6439	2.069	-1.761	0.078
Florida
regionFlorida Keys	0.0861	1.879	0.046	0.963
treatyi.i:regionFlorida Keys	-0.9606	3.184	-0.302	0.763
treatyi.w:regionFlorida Keys	5.1661	3.061	1.688	0.091
Variance Components
QE	504.3862
p value	0.0000
QM	9.3736
p value	0.0951

#save_kable(table3, "Table3_RegionalMod_16Jan20.pdf")

Supplementary Table S4. Temperature and region linear mixed effects model section using AICc. All models were run with random effects for study and species and a weight of 1/SE. The asterisk (*) denotes the selected model run for the final analysis.

warm.mag.tab <- AICc(warm.rate.region.mod, warm.rate.region2.mod, warm.rate2.mod, warm.rate.mod, warm.region.mod)

tableS4 <- kable(warm.mag.tab, digits = 2, booktabs = T) %>%
  kable_styling(font_size = 14, full_width = F)  %>% 
  row_spec(2, bold = T)
tableS4

	df	AICc
warm.rate.region.mod	7	109.99
warm.rate.region2.mod	8	75.39
warm.rate2.mod	7	72.85
warm.rate.mod	6	107.84
warm.region.mod	6	125.18

Supplementary Table S5. Temperature and region best fit linear mixed effects model output and 95% confidence intervals of the calcification rates in response to treatment temperature plotted in Figure 4.

tab_model(warm.rate.region2.mod, show.r2 = FALSE, show.aicc = TRUE)

	rate
Predictors	Estimates	CI	p
(Intercept)	0.74	-1.72 – 3.21	0.555
scale.temp	-0.23	-0.31 – -0.16	<0.001
region: Florida Keys	0.24	-1.05 – 1.53	0.713
I(scale.temp^2)	-0.36	-0.44 – -0.27	<0.001
days	0.01	-0.02 – 0.04	0.659
Random Effects
σ²	0.28
τ₀₀ _spec	0.20
τ₀₀ _study	0.18
ICC	0.57
N _study	5
N _spec	11
Observations	60
AIC	75.386

Supplementary Table S6. Aragonite saturation state and region linear mixed effects model section using AICc. All models were run with random effects for study and species and a weight of 1/SE. The asterisk (*) denotes the bet fit model.

arag.mag.tab <- AICc(acid.rate.region.mod, acid.rate.region2.mod, acid.rate2.mod, acid.rate.mod, acid.region.mod)

tableS6 <- kable(arag.mag.tab, digits = 2, booktabs = T) %>%
  kable_styling(font_size = 14, full_width = F)  %>% 
  row_spec(2, bold = T)
tableS6

	df	AICc
acid.rate.region.mod	7	57.97
acid.rate.region2.mod	8	49.85
acid.rate2.mod	7	47.98
acid.rate.mod	6	56.39
acid.region.mod	6	105.28

Supplementary Table S7. Aragonite saturation state (Ω_Arag) and region best fit linear mixed effects model output and 95% confidence intervals of the calcification rates in response to treatment Ω_Arag plotted in Figure 5.

tab_model(acid.rate.region2.mod, show.r2 = FALSE, show.aicc = TRUE)

	rate
Predictors	Estimates	CI	p
(Intercept)	-0.19	-3.06 – 2.69	0.900
scale.arag	0.10	0.07 – 0.14	<0.001
region: Florida Keys	0.47	-1.00 – 1.95	0.532
I(scale.arag^2)	-0.04	-0.07 – -0.02	0.001
days	0.02	-0.02 – 0.05	0.376
Random Effects
σ²	0.20
τ₀₀ _spec	0.17
τ₀₀ _study	0.25
ICC	0.68
N _study	5
N _spec	12
Observations	140
AIC	49.853

Supplemental Figures

oa.for <- forest(oa, slab=paste(OA$cite, OA$spec, OA$site, sep=" - "), order=order(OA$cite))

Supplementary Figure S1. Forest plot depicting individual standard mean difference (SMD) and confidence interval of each species-site-study combination included in the meta-analysis of ocean acidification only studies. Studies with confidence intervals not overlapping zero (dotted line) denote either a significant increase (positive values) or decrease (negative values) in calcification response under experimental acidification compared to the corresponding control treatment.

oa.fun <- funnel(oa, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)

Supplementary Figure S2. Funnel plot of the standard mean difference (SMD) against standard error of each species-site-study combination included in the meta-analysis of ocean acidification only studies. Dotted lines represent confidence intervals of all studies, background plot colour represents statistical significance of individual study, the black vertical line denotes the overall effect of all included ocean acidification studies.

pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS1_OAforest.pdf", width = 11, height = 12)
forest(oa, slab=paste(OA$cite, OA$spec, OA$site, sep=" - "), order=order(OA$cite))
dev.off()

pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS2_OAfunnel.pdf", width = 7, height = 6)
funnel(oa, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)
dev.off()

ow.for <- forest(ow, slab=paste(OW$cite, OW$spec, OW$site, sep=" - "), order=order(OW$cite))

Supplementary Figure S3. Forest plot depicting individual standard mean difference (SMD) and confidence interval of each species-site-study combination included in the meta-analysis of ocean warming only studies. Studies with confidence intervals not overlapping zero (dotted line) denote either a significant increase (positive values) or decrease (negative values) in calcification response under experimental warming compared to the corresponding control treatment.

ow.fun <- funnel(ow, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)

Supplementary Figure S4. Funnel plot of the standard mean difference (SMD) against standard error of each species-site-study combination included in the meta-analysis of ocean warming only studies. Dotted lines represent confidence intervals of all studies, background plot colour represents statistical significance of individual study, the black vertical line denotes the overall effect of all included ocean acidification studies.

pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS3_OWforest.pdf", width = 9, height = 10)
forest(ow, slab=paste(OW$cite, OW$spec, OW$site, sep=" - "), order=order(OW$cite))
dev.off()

pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS4_OWfunnel.pdf", width = 7, height = 5)
funnel(ow, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend="topleft")
dev.off()

com.for <- forest(both, slab=paste(IN$cite, IN$spec, IN$site, sep=" - "), order=order(IN$cite))

Supplementary Figure S5. Forest plot depicting individual standard mean difference (SMD) and confidence interval of each species-site-study combination included in the meta-analysis of the combination of ocean acidification and warming studies.

com.fun <- funnel(both, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)

Supplementary Figure S6. Funnel plot of the standard mean difference (SMD) against standard error of each species-site-study combination included in the meta-analysis of combined ocean acidification and warming only. Dotted lines represent confidence intervals of all studies, background plot colour represents statistical significance of individual study, the black vertical line denotes the overall effect of all included ocean acidification studies.

pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS5_bothforest.pdf", width = 9, height = 8)
forest(both, slab=paste(IN$cite, IN$spec, IN$site, sep=" - "), order=order(IN$cite))
dev.off()

pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS6_bothfunnel.pdf", width = 7, height = 6)
funnel(both, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)
dev.off()

sessionInfo()

## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] cowplot_1.0.0      kableExtra_1.1.0   sjlabelled_1.1.1  
##  [4] sjmisc_2.8.2       sjPlot_2.7.2       plotly_4.9.0      
##  [7] shiny_1.4.0        Rmisc_1.5          lattice_0.20-38   
## [10] MuMIn_1.43.6       RColorBrewer_1.1-2 tidystats_0.4     
## [13] forcats_0.4.0      stringr_1.4.0      dplyr_0.8.3       
## [16] purrr_0.3.3        readr_1.3.1        tibble_2.1.3      
## [19] tidyverse_1.2.1    lme4_1.1-21        plyr_1.8.4        
## [22] tidyr_1.0.0        ggplot2_3.2.1      MAd_0.8-2.1       
## [25] metafor_2.1-0      Matrix_1.2-17     
## 
## loaded via a namespace (and not attached):
##  [1] TH.data_1.0-10    minqa_1.2.4       colorspace_1.4-1 
##  [4] ellipsis_0.3.0    estimability_1.3  parameters_0.2.0 
##  [7] rstudioapi_0.10   ggrepel_0.8.1     mvtnorm_1.0-11   
## [10] lubridate_1.7.4   xml2_1.2.2        codetools_0.2-16 
## [13] splines_3.5.2     mnormt_1.5-5      knitr_1.25       
## [16] zeallot_0.1.0     jsonlite_1.6      nloptr_1.2.1     
## [19] ggeffects_0.12.0  Cairo_1.5-10      broom_0.5.2      
## [22] compiler_3.5.2    httr_1.4.1        sjstats_0.17.6   
## [25] emmeans_1.4.2     backports_1.1.5   assertthat_0.2.1 
## [28] fastmap_1.0.1     lazyeval_0.2.2    cli_1.1.0        
## [31] later_1.0.0       formatR_1.7       htmltools_0.4.0  
## [34] tools_3.5.2       coda_0.19-3       gtable_0.3.0     
## [37] glue_1.3.1        Rcpp_1.0.2        cellranger_1.1.0 
## [40] vctrs_0.2.0       nlme_3.1-141      crosstalk_1.0.0  
## [43] psych_1.8.12      insight_0.6.0     xfun_0.10        
## [46] rvest_0.3.4       mime_0.7          lifecycle_0.1.0  
## [49] MASS_7.3-51.4     zoo_1.8-6         scales_1.0.0     
## [52] hms_0.5.1         promises_1.1.0    parallel_3.5.2   
## [55] sandwich_2.5-1    yaml_2.2.0        stringi_1.4.3    
## [58] highr_0.8         bayestestR_0.4.0  boot_1.3-23      
## [61] rlang_0.4.1       pkgconfig_2.0.3   evaluate_0.14    
## [64] htmlwidgets_1.5.1 labeling_0.3      tidyselect_0.2.5 
## [67] magrittr_1.5      R6_2.4.0          generics_0.0.2   
## [70] multcomp_1.4-10   pillar_1.4.2      haven_2.1.1      
## [73] foreign_0.8-72    withr_2.1.2       survival_2.44-1.1
## [76] performance_0.4.0 modelr_0.1.5      crayon_1.3.4     
## [79] rmarkdown_1.16    grid_3.5.2        readxl_1.3.1     
## [82] data.table_1.12.6 digest_0.6.22     webshot_0.5.1    
## [85] xtable_1.8-4      httpuv_1.5.2      stats4_3.5.2     
## [88] munsell_0.5.0     viridisLite_0.3.0

Meta-analysis reveals reduced coral calcification under projected ocean warming but not under acidification across the Caribbean Sea

Colleen Bove, James Umbanhowar, and Karl Castillo